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High  Order  Accuracy  Computational  Methods 
In  Aerodynamics 
Using  Parallel  Architectures 
AFOSR  95-1-0074 

David  Gottlieb,  Chi-Wang  Shu,  P.f.  Fischer  and  1F.5.  Don 
Division  of  Applied  Mathematics,  Brown  University 

The  main  theme  of  this  research  is  the  application  of  high  order  accurate  schemes  to  complicated 
flow  problems.  The  advantage  of  using  high  order  schemes  for  long  time  simulations  is  widely 
recognized  by  now.  Also  for  problems  where  fine  details  of  the  flow  field  have  to  be  captured 
accurately  the  use  of  high  accuracy  schemes  is  mandatory.  These  two  classes  of  problems  encompass 
many  of  the  current  problems  in  scientific  computing. 

1.  Spectral  shock  capturing  techniques 

The  problem  in  applying  spectral  methods  to  shock  wave  problems  is  their  sensitivity  to  dis¬ 
continuities.  The  presence  of  shock  wave  creates  theoretical  as  well  as  practical  difficulties. 

Two  major  theoretical  breakthroughs  that  occured  in  the  last  period,  may  clear  the  way  to  the 
successful  implementation  of  spectral  methods  in  simulating  complicated  shock  waves  interactions. 

The  first  development  involved  the  solution  of  the  Gibbs  phenomenon.  The  Gibbs  phenomenon 
is  related  to  the  well  known  fact  that  the  rate  of  convergence  of  the  sequence  of  partial  sums  of  the 
Fourier  series  representation  of  a  function  deteriorates  rapidly  when  the  function  has  discontinuities. 
In  essence  the  new  result  says  that  the  first  N  Fourier  coefficients  of  a  piecewise  analytic  function 
(a  finite  number  of  bounded  jump  discontinuities  is  permitted)  carries  much  more  information  than 
can  be  revealed  by  forming  the  N  th  order  Fourier  partial  sum.  In  fact  it  has  been  shown  that  these 
N  coefficients  contain  enough  information  so  that  one  can  constructively  form  an  approximation 
to  the  unknown  piecewise  analytic  function  which  is  exponentially  accurate  in  N.  This  result  had 
been  extended  too  any  polynomial  method,  such  as  the  commonly  used  Chebyshev  method. 

The  second  development  concerns  the  convergence  of  spectral  methods  for  non-linear  hyperbolic 
equations.  It  has  been  shown  that  with  a  suitable  addition  of  (spectrally  amaU)  artificial  dissipation 
the  method  converges.  Moreover  it  has  been  shown  that  the  usual  filtering  technique  can  be  cast 
in  terms  of  the  above  analized  artficial  dissipation. 

While  the  above  results  create  a  breakthrough  that  indicate  that  high  order  methods  are  useful 
for  shock  wave  calculations,  a  lot  of  work  has  still  to  be  .  done  in  making  the  process  efficient. 
Stabilizing  the  scheme  in  an  optimal  way,  and  extracting  the  information  in  an  efficient  way  have 
to  be  further  studied. 


2.  High  order  ENO  schemes 

We  have  been  continuing  our  investigation  of  ENO  schemes  using  point-flux  and  TVD  Runge- 
Kutta  time  discretization  formulations  in  the  following  directions: 

(1)  Towards  the  convergence  issues,  we  have  investigated  entropy  consistency  for  high  order 
Hermite  type  finite  difference  and  discontinuous  Galerkin  methods  as  a  first  step.  We  have  been 
able  to  obtain  a  cell  entropy  inequality  for  the  square  entropy,  for  such  high  order  schemes  (no 
restriction  on  order  of  accuracy)  without  resorting  to  the  help  of  nonlinear  limiters.  The  result  is 
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valid  for  multi  space  dimensions  with  arbitrary  triangulations  and  for  any  fluxes  (no  restriction  on 
convexity).  This  is  a  significant  improvement  in  obtaining  cell  entropy  inequalities  since  all  the 
previous  work  in  this  direction  must  either  resort  to  modifications  to  existing  limiters  or  resort  to 
complicated  global  analysis,  and  be  restricted  to  second  order  accuracy.  We  plan  to  continue  the 
investigation  towards  full  convergence  proofs; 

(2)  Application  of  ENO  schemes  to  combustion  problems.  The  first  test  case  is  shock  interaction 
with  hydrogen  bubbles.  Different  configuration  of  bubbles  is  studied  to  see  the  effect  of  shock 
interaction.  High  order  accuracy  is  crucial  in  this  problem  due  to  the  detailed  structures  of  the 
solution  behind  the  shock.  This  project  is  on  going.  Stiff  source  terms  must  be  treated  adequately 
for  future  tests; 

(3)  Comparison  of  two  different  formulation  of  ENO  schemes:  finite  volume  vs  finite  difference. 
The  comparison  is  performed  on  problems  related  to  curved  boundary,  inflow  outflow  boundary 
conditions,  shocks  oblique  to  the  grids,  and  CPU  time.  It  is  found  out  that  finite  difference  version 
of  ENO  scheme  (the  one  based  on  point  values  and  numerical  fluxes  of  Shu  and  Osher)  is  much 
cheaper  to  run,  and  can  obtain  results  comparable  to  finite  volume  ENO  in  most  of  the  test  cases. 

(4)  ENO  schemes  based  on  rational  function  building  blocks  are  being  investigated.  As  a  first 
step  we  are  investigating  approximation  results.  The  potential  here  is  that  in  most  rapid  transition 
regions,  rational  functions  approximate  the  true  solution  better  than  polynomials.  This  project  is 
on  going. 

The  results  of  comparing  the  two  formulations  of  ENO  schemes  (the  cell- averaged  version  and 
the  point  value  version  indicate  that  for  most  test  cases,  the  two  formulations  of  ENO  schemes 
yields  the  same  accuracy  whereas  the  point  value  ENO  scheme  is  much  faster. 

3.  Parallel  computing  for  high  order  schemes 

(a)  Spectral  Elements  for  incompressible  flows 

The  research  objective  here  is  to  develop  parallel  algorithms  for  computational  fluid  dynam¬ 
ics  (CFD)  which  will  permit  solution  of  incompressible  flows  with  the  accuracy  and  resolution 
demanded  by  large  eddy  simulations  (LES)  of  turbulent  flows  in  complex  geometries.  The  work 
is  subdivided  into  three  principal  areas:  high-order  flexible  discretizations,  fast  multi-level  itera¬ 
tive  solvers,  and  implementation  of  LES  modeling  technology.  All  of  this  work  is  being  developed 
within  the  computational  framework  of  distributed  memory  architectures  which  provide  a  favorable 
price/performance  ratio  for  this  class  of  problems. 

Our  numerical  approach  is  based  upon  the  spectral  element  method,  which  retains  many  of 
the  essential  features  of  global  spectral  discretizations,  namely,  rapid  (exponential)  convergence, 
minimal  numerical  dissipation  and  dispersion  per  degree-of-freedom,  and  efficient  tensor  product 
factorization  of  spatial  derivative  operators.  The  computational  domain  is  subdivided  into  large 
macro-elements  and  the  solution,  data,  and  geometry  within  each  element  are  expressed  in  terms 
of  high-order  polynomials.  The  use  of  a  weighted- residual  procedure  permits  a  reduction  in  inter¬ 
element  continuity  requirements  from  to  C°,  which  in  turn  leads  to  a  reduction  in  inter-processor 
communication.  The  locally-structured/globally- unstructured,  approach  of  the  spectral  element 
discretization  is  ideally  suited  to  the  two-level  memory  hierarchy  associated  with  distributed  mem¬ 
ory  parallel  computers. 

The  performance  of  general  geometry  incompressible  codes  is  largely  tied  to  the  speed  of  the 
elliptic  solver  for  the  pressure;  development  of  fast  multi-level  iterative  solvers  is  a  major  focus 
of  our  current  efforts.  Presently,  we  employ  a  two-level  conjugate  gradient  based  solver  which 
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uses  deflation  to  project  out  the  coarse  grid  modes,  thereby  reducing  the  condition  number  of  the 
underlying  iteration  matrix.  The  system  is  preconditioned  by  local  finite-element  based  operators 
which  are  significantly  less  expensive  to  invert  than  their  high-order  counterparts.  In  effect,  one 
can  obtain  high-order  accuracy  at  low-order  cost.  We  have  addressed  all  of  the  parallel  issues 
associated  with  this  approach,  including  the  use  of  an  efficient  parallel  coarse-grid  solve,  and  are 
able  to  compute  Navier-Stokes  solutions  at  a  resolution  of  3  million  gridpoints  on  the  512  node 
Intel  Delta  at  the  rate  of  4  minutes  per  time  step  (6.5  GFLOPS).  To  improve  upon  this  result, 
we  are  currently  developing  a  preconditioner  based  upon  a  Schur-complement  formulation  for  the 
interface  variables.  This  approach  leads  not  only  to  better  conditioning,  but  also  to  reduced  degrees- 
of-freedom,  which  in  turn  permits  the  use  of  projection  techniques  in  which  very  good  initial  guesses 
can  be  generated  by  projecting  the  residual  onto  results  from  previous  time  steps. 

Accurate  numerical  simulation  of  many  high  Reynolds  number  engineering  flows  will  continue 
to  be  limited  by  resolution  requirements  for  the  foreseeable  future.  However,  progress  is  being  made 
to  the  point  where  the  combination  of  high-resolution  (  >  10®  gridpoints  )  and  advanced  large-eddy 
simulation  (LES)  models  will  be  able  to  conquer  many  important  problems  in  the  near  future.  Our 
goal  is  to  couple  the  latest  LES  technology  developed  within  the  turbulence  community  with  a 
general  geometry  Navier-Stokes  code  capable  of  solving  these  demanding  problems. 

Because  of  the  large  number  of  degrees-of-freedom  involved,  effective  use  of  high-performance 
distributed  memory  parallel  architectures  is  essential  to  economic  resolution  of  these  problems. 
Principal  parallel  issues  to  be  addressed  include:  development  of  coarse-grid  solve  strategies  which 
will  remain  competitive  as  the  number  of  processors  and  dimension  of  the  coarse-grid  system  con¬ 
tinue  to  increase;  and  development  of  optimal  communication  strategies  for  the  complex  subdomain 
interfaces  arising  from  nonconforming  discretizations. 


(b)  Spectral  methods  for  compressible  flows. 


Spectral  methods  involve  the  approximation  of  the  unknown  solution  in  terms  of  global  poly¬ 
nomials.  This  fact  make  them  difficult  to  implement  on  parallel  computers.  A  popular  method  to 
overcome  the  limitations  of  spectral  methods  is  to  use  multidomain  techniques,  in  which  a  complex 
domain  is  decomposed  into  several  simpler  subdomains.  This  method  has  been  applied  successfully 
to  incompressible  flows  (the  Specral-Element  technique)  or  to  problems  in  structural  mechanics 
(the  h-p  method). 

Multidomain  spectral  methods  are  suitable  for  coarse  grain  parallel  computing,  each  domain  is 
assigned  to  a  different  processor. 


The  main  question  is:  Are  multidomain  methods  efficient?  This  question  hcis  not  been  yet 
answered  for  those  methods  applied  to  hyperbolic  equations.  If  we  denote  by  W{p^N)  the  work 
involved  in  approximating  k  waves  using  p  sub  domains  (and  N  points  in  each  domain)  to  obtain 
an  error  of  at  most  e~'^,  then  W(p,  N)  is  minimized  if 
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Thus  the  optimal  number  of  subdomain  increases  with  the  complexity  of  the  problem  (or  number  of 
waves)  but  decreases  if  the  required  accuracy  increases.  This  result  may  serve  as  a  guideline  to  the 
optimal  number  of  subdomain.  The  formula  above  can  be  suitably  modified  for  parallel  computers. 


A  key  issue  in  the  application  of  multidomain  spectral  methods  to  the  numerical  solution  of 
hyperbolic  equations  is  the  interface  boundary  conditions.  This  leads  to  the  question  of  the  imposi¬ 
tion  of  boundary  conditions,  both  analytic  and  numerical,  in  the  numerical  solutions  of  systems  of 
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hyperbolic  equations.  For  truly  time  dependent  problems  stability  in  the  classical  sense  (Lax  and 
G-K-S  stability)  is  not  enough.  Even  stable  schemes  may  exhibit  a  non-physical  growth  in  time. 
From  a  practical  point  of  view,  in  order  to  achieve  reasonable  accuracy  for  large  time,  meshes  much 
too  fine  for  the  computers  available  in  the  foreseeable  future  are  required. 

We  have  found  a  systematic  procedure  for  designing  time-stable,  as  well  as  G-K-S  stable  schemes 
of  high-order  accuracy.  The  new  schemes  are  guaranteed  to  be  time-stable  for  any  hyperbolic  system 
(as  long  as  the  system  has  a  bounded  energy).  We  have  extended  this  methodology  to  Navier  Stokes 
equations  in  three  space  dimensions.  We  have  showed  that  the  SAT  boundary  condition  assures 
the  correct  behavior  as  the  Reynolds  number  tends  to  infinity. 

(c)  Parallel  implementation  of  ENO  schemes  on  CM-5 

The  main  cost  in  ENO  schemes  is  in  its  logic  step  in  choosing  local  stencils  by  comparing  divided 
difference  tables  of  the  function.  Although  great  effort  has  been  spent  on  efficiently  vectorizing  this 
part  for  CRAY  supercomputers,  due  to  the  inevitable  gathering-scattering  process,  ENO  schemes 
still  do  not  run  very  fast  on  CRAY  computers.  Recently  we  have  been  exploring  the  structure 
of  the  ENO  algorithm  to  suit  the  parallel  structure  of  CM-5.  The  algorithm  has  been  slightly 
re-formulated  (in  a  mathematically  equivalent  way)  to  reduce  communications  and  to  eliminate 
communications  between  other  than  next  neighbors,  at  the  price  of  a  slightly  increased  operation 
count.  Our  CM-5  two  dimensional  ENO  code  for  compressible  Euler  and  Navier-Stokes  equations 
is  a  magnitude  faster  than  our  ENO  code  on  CRAY  for  a  400x400  grid.  Three  dimensional  code 
also  shows  a  speed  up,  although  it  has  not  been  optimized  yet  (for  three  dimensions,  storage  is  a 
big  restriction,  and  the  structure  of  the  program  must  be  modified  accordingly).  Currently  we  are 
trying  to  improve  the  CM-5  code  and  applying  it  to  reactive  flows.  . 

4.  Fuel  air  mixing  enhanced  by  shock  induced  vortices 

In  designing  supersonic  combuster  for  the  next  generation  of  supersonic  transport,  we  look  for 
an  efficient  combuster  such  that  : 

•  allow  better  load  to  weight  ratio  by  carrying  less  fuel, 

•  reduce  chemical  product  that  contribute  to  pollution  due  to  incomplete  burning  of  fuel. 

One  of  the  technique  currently  under  study  in  enhancing  mixing  of  a  hydrogen  jet  (fuel)  and 
air  (oxidizer)  is  to  allow  the  existing  shock  inside  the  combustion  chamber  to  interact  with  the 
hydrogen  jet  .  By  doing  so,  vorticity  are  generated  according  to  the  vorticity  equation  of  the 
Navier  Stokes  equation.  The  pressure  gradient  of  the  shock  and  the  density  gradient  between  the 
air  and  hydrogen  provides  an  efficient  mechanism  for  vorticity  generation  along  the  surface  of  the 
hydrogen  jet.  The  vorticity  forced  the  jet  to  curl  up  and  stressed.  The  increased  surface  area  allows 
greater  mixing  of  air  and  hydrogen  where  combustion  can  take  place.  Even  though  the  problem 
is  three  dimensional  steady  state,  it  had  been  argued  that  the  three  dimensional  steady  flow  can 
be  directly  associated  with  a  corresponding  two  dimensional  unsteady  calculation  in  this  particular 
physical  setting. 

Numerical  simulations  are  a  necessary  part  of  this  investigation.  For  these  type  of  problem,  it 
is  important  to  capture  the  complex  physics  with  high  accuracy.  Finite  difference  scheme  ,  with 
their  inherent  dissipative  nature  for  stability  reason,  can  only  yields  a  quantitative  result  of  the 
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flow  fields.  Often,  for  long  time  integration,  the  loss  of  accuracy  in  the  earlier  stage  affect  the 
development  of  the  flow  in  later  time. 

We  had  applied  two  codes  for  simulating  the  above  problem.  A  spectral  shock  capturing  code 
and  a  finite  difference  ENO  code.  We  have  employed  all  the  theoretical  techniques  developed  in 
order  to  run  succesfuUy  the  codes. 

we  computed  the  solution  of  these  problem  with  various  configuration  of  multiple  jets  place¬ 
ment.  Our  results  found  that  different  initial  placement  of  the  jets  yields  distinct  final  configurations 
(Moreover,  the  spectral  code  is  fully  capable  of  capturing  the  fine  scale  structure  of  the  interac¬ 
tions,  including  some  that  had  not  been  seen  in  finite  difference  code.  One  distinct  feature  of  this 
calculation  is  the  penetration  of  an  air  stream  (heavy  fluid)  into  the  hydrogen  (light  fluid)  causing 
instability  that  exhibit  a  mushroom  shape  structure  This  study  can  be  used  to  guided  researcher 
in  develop  fuel  jet  configuration  for  the  scramjet  engine  with  greater  confident.  Our  goal  in  this 
research  is  to  extend  this  problem  to  a  three  dimensional  steady  and/or  unsteady  simulation  with 
full  chemistry  model. 


5.  Multiscale  Computing 

Many  relevant  physical  phenomena  involve  infinite  number  of  scales.  In  nonlinear  problems  it 
is  desirable  to  find  a  way  to  take  into  account  the  effect  of  the  scales  which  are  neglected  on  those 
which  are  taJcen  into  account. 

Two  approaches  have  been  investigated  for  this  type  of  problems:  Nonlinear  Galerkin  Methods 
and  Wavelet  based  Schemes. 

We  have  studied  the  implementation  of  the  Nonlinear  Galerkin  in  the  context  of  collocation 
discretizations.  We  have  found  interesting  characterizations,  in  the  physical  space,  of  a  small  scale 
function  and  a  large  scale  function. 

Using  the  above  we  have  presented  an  efficient  pseudospecral  NLG  scheme  for  the  periodic 
Burgers  equation.  The  case  of  Chebyshev  approximations  to  nonperiodic  problems,  in  which  the 
concept  of  large  and  small  scales  have  to  be  redefined  is  currently  under  investigation. 

In  a  thesis  by  L.  Jameson  and  in  a  subsequent  paper  ,  the  Daubechies  wavelet  based  differenti¬ 
ation  matrix  was  constructed.  The  relationship  between  this  matrix  and  finite  difference  methods 
waa  clarified.  This  serves  as  a  basis  for  current  work  by  Bruce  Bauer  a  doctoral  student  on  wavelet 
optimized  finite  difference  methods. 

Jameson  also  analyzed  the  differentiation  matrix  based  on  the  compactly  supported  Daubechies 
wavelets.  He  showed  that  in  this  case  there  is  a  loss  of  the  superconvergence.  The  same  result  holds 
for  the  finite  element  schemes. 
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